#include<iostream>
#include<cmath>
#include<vector>
#include "../programming1/EquationSolver.h"
#include "LUFactorization.h"
using namespace std;

class CubicBSpline{
private:
	int L, R;
	vector <double> coefficient;
public:
	CubicBSpline(const int& x0, const int& n, const vector<double>& f, const string& mode = "natural", const double& m_0 = 0, const double& m_n = 0) : L(x0) {
		R = L + n;
		vector <double> d(n+1), u(n), l(n), b(n+1);
		for (int i = 1; i < n; ++ i) {
			d[i] = 4; //对角
			l[i-1] = 1; //上面
			u[i] = 1;//下面
			b[i] = 6 * f[i];
		}
		if (mode == "natural" || mode == "complete") {
			d[0] = 2, u[0] = 1, b[0] = 3 * f[0] + m_0;
			d[n] = 2, l[n-1] = 1, b[n] = 3 * f[n] - m_n;
			vector <double> t = LU(d, u, l, b);
			coefficient.resize(n+3);
			for (int i = 0; i <= n; ++ i) coefficient[i+1] = t[i];
			coefficient[0] = coefficient[2] - 2 * m_0;
			coefficient[n+2] = coefficient[n] + 2 * m_n;
		}
		else if (mode == "ssd") {
			d[0] = 6, b[0] = 6 * f[0] - m_0;
			d[n] = 6, b[n] = 6 * f[n] - m_n;
			vector <double> t = LU(d, u, l, b);
			coefficient.resize(n+3);
			for (int i = 0; i <= n; ++ i) coefficient[i+1] = t[i];
			coefficient[0] = m_0 - coefficient[2] + 2 * coefficient[1];
			coefficient[n+2] = m_n - coefficient[n] + 2 * coefficient[n-1];
		}
	}
	double getValue(const double& _x)  {
		if (_x < L || _x > R) throw "Out of Range!";
		int i = floor(_x);
		double res = 0;
		if(i-2 >= L-2 && i-2 <= R) res += coefficient[i-2-L+2] * B(3, i-2, _x);
		if(i-1 >= L-2 && i-1 <= R) res += coefficient[i-1-L+2] * B(3, i-1, _x);
		if(i   >= L-2 && i   <= R) res += coefficient[i  -L+2] * B(3, i  , _x);
		if(i+1 >= L-2 && i+1 <= R) res += coefficient[i+1-L+2] * B(3, i+1, _x);
		return res;
	}
};


CubicBSpline CubicBSplineInterpolation(Function & f, const int& l, const int& n, const string& mode = "natural") {
	vector <double> a(n+1);
	for(int i = 0; i <= n; ++ i) a[i] = f(l+i);
	double m_0, m_n;
	if (mode == "natural") m_0 = m_n = 0;
	else if (mode == "complete") m_0 = f.diff(l), m_n = f.diff(l+n);
	else if (mode == "ssd") m_0 = f.diff2(l), m_n = f.diff2(l+n);
	return CubicBSpline(l, n, a, mode, m_0, m_n);
}
